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Abstract 

Motivated by recent advances in vesicle engineering, we consider theoretically the locomotion of shape- 
changing bilayer vesicles at low Reynolds number. By modulating their volume and membrane composition, 
the vesicles can be made to change shape quasi-statically in thermal equilibrium. When the control param- 
eters are tuned appropriately to yield periodic shape changes which are not time-reversible, the result is a 
net swimming motion over one cycle of shape deformation. For two classical vesicle models (spontaneous 
curvature and bilayer coupling), we determine numerically the sequence of vesicle shapes through an en- 
thalpy minimization, as well as the fluid-body interactions by solving a boundary integral formulation of the 
Stokes equations. For both models, net locomotion can be obtained either by continuously modulating fore- 
aft asymmetric vesicle shapes, or by crossing a continuous shape-transition region and alternating between 
fore-aft asymmetric and fore-aft symmetric shapes. The obtained hydrodynamic efficiencies are similar to 
that of other low Reynolds number biological swimmers, and suggest that shape-changing vesicles might 
provide an alternative to flagclla-based synthetic microswimmcrs. 
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I. INTRODUCTION 



The preeminence of viscous dissipation over inertial effects at low Reynolds numbers leads to 
many interesting consequences for life and engineering efforts at the micron-scale. In particular, 
swimming at zero Reynolds number is impossible using time-reversible motions, a result known as 
the Scallop theorem [1]. As a result, at least two actuation degrees of freedom are necessary to 
generate locomotion. The breaking of this time-reversal symmetry has been studied both from a 
mathematical point of view, and in the context of modeling real organisms [2-6]. Unlike in high 
Reynolds number flows, such as those relevant in describing the swimming of fish and flying of 
birds, fluid motion at low Reynolds numbers is set almost instantaneously by the time-dependent 
geometries of the immersed bodies. Thus it is natural to inquire about the shapes of immersed 
(and possibly fluctuating) cell membranes, and their relationships to locomotion. 

Membranes composed of lipid bilayers are ubiquitous in nature, and the study of bilayer vesicles 
as a model system for biological cells has yielded significant insight into their behavior [7, 8]. 
In addition to the biological relevance of lipid bilayer vesicles, or liposomes, advances in self- 
assembly have paved the way for other types of vesicles to be developed experimentally [9, 10]. 
Vesicles assembled from block copolymers [11], liquid crystal amphiphiles [12], and membranes 
with embedded proteins or anchored polymers [13-17] all have tunable material properties which 
can be manipulated with unprecedented control [18, 19]. It is also well known that many biological 
cells actively modify or maintain the shapes of their membranes [20, 21], either for developmental 
[22] or locomotive processes [23, 24]. 

Recently, synthetic microswimmers inspired by the locomotion of eukaryotic cells have been 
successfully designed in experiments [25], exploiting the planar beating of a flagellum-like organelle. 
Beyond biomimetic engineering, other small-scale synthetic swimmers or swimming strategies have 
also been proposed, both theoretically and experimentally [1, 6, 26-32]. One recently-studied 
example is a self-propelled colloidal particle which exploits asymmetrically-distributed chemical 
reactions to swim in a viscous fluid [33, 34]. 

In the same spirit, we consider theoretically in this paper a novel swimming mechanism based on 
prescribed shape transformations of a bilayer vesicle. By modulating only its volume and membrane 
composition, the vesicle can be made to change shape quasi-statically in thermal equilibrium. For 
two different theoretical vesicle models, we determine numerically the vesicle shapes through an 
enthalpy minimization, and the fluid-body interactions by solving a boundary integral formulation 
of the Stokes equations. When the control parameters are tuned appropriately to yield periodic but 
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not time-reversible shape changes, we show that net locomotion can be obtained. Swimming arises 
either by continuously modulating fore-aft asymmetric vesicle shapes, or by crossing a continuous 
shape-transition region and alternating between fore-aft asymmetric and fore-aft symmetric shapes. 
In addition, the calculated hydrodynamic efficiencies are shown to be similar to that of other 
common low Reynolds number propulsive mechanisms. 

Our paper is organized as follows. We begin with a general discussion of the practical realiza- 
tion of controlled shape-changing vesicles, in particular the relevant time scales, and the possible 
actuation mechanisms. Two classical curvature-mediated vesicle models (spontaneous curvature 
and bilayer coupling) are presented, and the formulations used for the shape calculation and the 
numerical fluid-interaction model are introduced. We then discuss examples of vesicle shape cycles 
that yield a swimming motion, examine the fluid flow that develops around the vesicles during 
their deformation cycles, and compute the corresponding swimming speeds and hydrodynamic 
efficiencies. 

II. A ROADMAP TO VESICLE LOCOMOTION 

A vesicle immersed in a viscous fluid experiences a highly coupled array of forces, such as 
those generated by membrane tension, internal pressure, membrane (bending) elasticity, and the 
surrounding viscous fluid dynamics. In the microscopic, viscous environments relevant to our 
consideration, the Reynolds number, Re, is very small: Re = pU c L c / [i <C 1, where p is the fluid 
density, ji is the fluid shear viscosity, and U c and L c are characteristic velocity and length scales 
of the vesicle. The fluid behavior at low Reynolds number is highly dependent upon the immersed 
boundary geometry, and the resultant forces include not only local, but also non-local responses 
to its motion. 

A general study of vesicle dynamics should take non-equilibrium shapes into account, as even 
simple liposomes that can be created in situ can interact relatively quickly with the environment. 
It is possible to design experiments where carefully constructed initial conditions and lipid species 
lead to equilibrated vesicle shapes that are non-trivial, but in order to apply morphological changes 
and induce locomotion, a reversible parameter-changing mechanism is desirable. 

For our first approach to vesicle swimming, we consider in this paper a "stiff membrane" regime. 
The characteristic time of membrane relaxation in a viscous fluid is given by t re i = fir^/n, where k 
is the elastic bending modulus of the membrane, and tq is a characteristic radius of curvature. If 
we choose the maximum radius of the vesicle for the characteristic length scale L c , then ro < L c . 
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For parameter variation significantly slower than the membrane relaxation rate, i.e. for a cycle 
time scale t cyc i e 3> t re i, then we operate safely within the decoupled regime. In this case, we may 
thus assume that there are no hydrodynamically induced shape changes, and that the shapes are 
determined quasi-statically in equilibrium. Using this time scale t re i, we can also set a maximum 
swimming velocity scale, U c = K/fir^. Similar scaling arguments have been made in Refs. [35, 36]. 
For biologically relevant systems in water, k ~ 100 kgT, \i ~ 10~ 3 Pa s, p = 1 g/cm 3 , and 
ro < 1 — 10 /im, leading to L c 1 — 10 fim, U c « 1 — 10 /im/s, t re \ « 0.01 — Is and Re ps 10 -4 . 
For a vesicle with length scale L c = 10 /im, diffusive time scales are approximately 10 4 s, and thus 
negligible for the time being. In addition, we neglect thermal fluctuations in the determination of 
the vesicle shape, as they come in as a perturbation about the mean equilibrium shape of order 
(k B T/K) 2 , which is very small under most conditions [7]. 

There are a number of different physical means by which a vesicle shape can be changed in a 
controlled fashion, and the methods could be different depending on the type of vesicle considered. 
We will consider two such means, internal volume changes and local membrane compositional 
changes. 

One experimentally feasible example of a possible volume-changing mechanism is a light-induced 
osmotic change. In an ordinary biological membrane the bilayer is embedded with numerous 
proteins, many of which are sensitive to mechanical forces, chemical gradients, or light. The 
protein bacteriorhodopsin, for example, is sensitive to green light, and in response to a signal the 
protein opens and closes like a valve [37]. The presence of such ion channels or active proteins on the 
surface of a membrane can cause osmotic changes of the fluid volume contained within the vesicle 
[19]. Recently, vesicle volume control was demonstrated via pH modulation of block copolymer 
networks along the surface of membrane [38] . The vesicles in this study were well separated from 
regimes associated with morphological transition, and thus changes in osmotic pressure induced 
only a volume change, leading to a "breathing" vesicle. 

Adjusting the membrane composition requires a more indirect experimental approach. Some 
bilayers are composed of different species of constituent parts, leading to an inherent mismatch 
between the intrinsic curvatures. In other words, there is an intrinsic curvature that would develop 
across the bilayer in the absence of other considerations. Because of the inherent difficulty in 
measuring these quantities it is likely to be more difficult to specify an exact change from one value 
of intrinsic curvature to another. However, the actual process of changing the intrinsic curvature 
can be achieved through inducing chemical changes of the lipid constituents of the membrane [38] , 
or by conformational changes of polymers grafted to the surface of the vesicle [39]. 
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FIG. 1: Schematic illustration of a possible control mechanism for vesicle shape-change and swimming: 
Axisymmetric bilayer vesicle with embedded reactive polymers and polymers grafted to its surface, (a) — > 
(b): Short frequency light impinges on the vesicle, catalyzing a de-polymerization reaction amidst the 
particle chains, and increasing the fluid volume available to the vesicle, (b) — > (c): A second frequency 
of light induces the grafted polymers to coil up, inducing an entropic repulsion from the membrane and 
changing the macroscopic morphology, (c) — > (d): The dispersed particles begin to polymerize back to 
their initial configuration, deflating the vesicle, (d) — > (a): A third frequency of light is used to uncoil the 
polymers, relaxing the entropically induced curvature and returning the vesicle to its initial state. 

By combining two shape-changing mechanisms, it would in theory be possible to achieve a 
periodic shape cycle which is not time-reversible, yielding a net locomotion. One of many possible 
configurations that could produce a cycle in shape space is displayed schematically in Fig. 1, where 
we consider a bilayer vesicle with embedded reactive polymers and with polymers grafted to its 
surface. In the first step (Fig. la— )-b), a photo-chemical polymerization reaction is catalyzed by 
green (short wavelength) light, and the polymer chains in the interior of the vesicle disperse into 
a solution of particles, thus increasing the available volume within the vesicle. At a later time, 
another frequency of light (red, or longer wavelength) impinges on the vesicle, and the grafted 
polymers change from a distended to a coiled conformation, inducing an entropic repulsion and 
changing the curvature of the membrane (Fig. lb— >c). Over time the dispersed particles will 
polymerize and return the vesicle to its original volume (Fig. lc— >d), and finally a third frequency 
of light (blue or very short wavelength) can be used to change the conformation of the polymers 
to distended once more, returning the vesicle to its original state (Fig. Id— >-a). 

While osmotic volume change or chemical-induced composition alteration are two possible ex- 
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perimental methods, not only these examples in no way constitute the full set of possibilities, but 
also they might be difficult to implement experimentally. Other experimental techniques already 
exist (see Refs. [19, 38]), or may be developed in the near future that could be more suited for 
controlled two-parameter change. 

Rather than suggest specific experimental methodology whose specifics would depend not only 
on the particular material of the bilayer vesicles, but also on the parameter alternation methods, 
we adopt in this paper a simplified modeling approach that highlights the qualitative pieces that 
are required in order to transform a motionless vesicle into a locomotive cargo-carrier. In parallel 
to the the various practical mechanisms that could be used to implement such shape changes 
experimentally, it is of fundamental interest to ask theoretically the question of prediction and 
performance. Would shape change indeed lead to locomotion of the vesicle? How efficient would 
it be? Can we quantitatively predict the resulting swimming speed and the work done against 
the fluid to achieve it? This is the approach taken in this paper. Considering two simplified 
vesicle models, and for slow modulations of the vesicle shapes, we introduce below a computational 
framework able to quantitatively predict swimming kinematics and performance. 

III. DYNAMICS OF COUPLED FLUID-BODY SYSTEM 
A. Vesicle physics 

While real biological membranes have multiple constituents, all interacting in non-trivial ways, 
minimal models can still help to illuminate the fundamental physics of such systems. For length 
scales on which a membrane is approximately flat a Monge parameterization can be employed [40— 
42], but for a closed bilayer vesicle the curvatures can become very large and the small geometric 
gradient assumption may break down. In order to characterize the shapes of such objects, an 
enthalpy must be extremized and the full nonlinear shape equations so generated must be solved. 
There are many models that could be used to describe the physics of curvature-mediated vesicle 
morphology. In this paper we will consider two classical models as case studies. These formulations, 
known respectively as the spontaneous curvature and bilayer coupling models, have both been 
used in classical work [43] and correspond to different interaction dynamics between the membrane 
monolayers. Both of these models also include exactly two free parameters, which enable us to 
explore the breaking of the Scallop theorem, and the generation of locomotion via a change in 
morphology. 
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The enthalpy functional, F, in the spontaneous curvature model takes the following form [43] 

F = ^l (C l + C 2 -C Q fdS + Y,A + PV, (1) 
1 JS(t) 

where C\ and C 2 are the principal membrane curvatures, and £ and P are Lagrange multipliers 
which constrain the surface area A and volume V (physically they correspond to the membrane 
tension and pressure difference across the interface). In Eq. (1), S(t) denotes the time-dependent 
surface boundary, and Cq is the spontaneous curvature, which introduces an inherent mismatch in 
equilibrium preference of the membrane curvature. This quantity along with a fixed volume and 
surface area completely specifies the ensemble. Thus the spontaneous curvature model has area, 
volume, and integrated spontaneous curvature constrained, and we select as the control parameters 
the volume V and the spontaneous curvature Co (the fixed surface area merely selects the overall 
size of the vesicle). 

In contrast, in the bilayer coupling model, the enthalpy functional G assumes the area difference 
AA between the membrane monolayers to be constant. One possible representation of this area 
difference is in terms of the integrated mean curvature, 

M= [ (C 1 + C 2 )dS. (2) 

Js(t) 

Then the area difference is AA = 2hM + 0(h 2 /A), where h is the distance between monolayers 
[43]. The enthalpy then takes the form 

G=l i! (Ci + C 2 ) 2 dS + Y! A + PV + Q M, (3) 
2 Js(t) 

where P, and Q are Lagrange multipliers associated with A (area), V (volume) and M (in- 
tegrated mean curvature) respectively. We select as control parameters the volume V and the 
integrated mean curvature M. 

It is important to note that the functionals F and C are related via a Legendre transform, 
(X', Q) — > (S + kCq/2, —2k Co), and thus describe the same system in a different ensemble. Phys- 
ically, the spontaneous curvature model corresponds to a bilayer in which the monolayer admits 
stretching or compression during bending, and thus finds an equilibrium distribution that has a 
preferred curvature. If the bilayer is composed of more than one species of lipid, each of which has 
a different preferred curvature (i.e. radius of gyration), it is likely that the membrane will actually 
prefer to be in a non-flat state. Conversely, the bilayer coupling model corresponds to a system 
that enforces that both monolayers are incompressible. The area difference between monolayers 
stays approximately constant on the timescales relevant to our consideration, and as long as the 
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FIG. 2: Parameterization of an axisymmetric bilayer vesicle. We assume axisymmetry about the z-axis. 
The surface is described by x = (r(s, t), z(s, t)) in cylindrical coordinates, with s an arc-length parameter, t 
the unit tangent vector, n the outward pointing normal vector, and i\) the angle between the x-axis and t. 

distance between layers remains very small this implies that the integrated mean curvature also 
remains constant. 



B. Determination of the vesicle shape 

Assuming an axisymmetric vesicle shape, the body surface S(t) is parameterized at each time t 
as illustrated in Fig. 2. The arc-length measured along the surface in the x-z plane is denoted by 
s € [0,1/], with t the unit tangent vector, n the outward pointing normal vector, and ip the angle 
between the x-axis and t. The body surface is represented in cylindrical polar coordinates, 

x(s, 4>, t) = x(s, (f>, t) + z (t)z = (r(s, t) cos((f>),r(s, t) sin(</>), z(s, t) + z (t)), (4) 

where <f> € [0, 2tt) is the azimuthal angle, the surface x is taken to have its center of volume at 
the origin, and zo(t) is a translation of that center of volume which depends upon the fluid in- 
teraction. Under this parameterization, the principal membrane curvatures are C\ = dip/ds and 
C2 = sinip/r. Upon insertion into either of the enthalpy functionals F or G, and performing a vari- 
ational extremization, we obtain the following system of first-order ordinary differential equations 
to describe the energetically stationary vesicle shapes at time t [43] 

^ = K, (5) 

K 1 , T • / 1 cos^sin^ 1 

K s = cos^H — smip -\ 7L 1 — Frcosip, (oj 

r r r 2 

(K-Cq) 2 sin 2 ^ . 
7s = ~ + Prsmip + T,, (7) 

r s = cos tp. (8) 

Here K is an auxiliary function used to make the system of equations first-order (physically it 
corresponds to the curvature), 7 is the Lagrange multiplier that enforces the interdependence of tp 



and r, and the subscript s denotes a derivative with respect to the arclength. The vesicle shape at 
time t is set by Eqs. (5-8), subject to the four boundary conditions r(0,t) = r(L,t) = ip(0,t) = 
and ip(L,t) = ir. Once the angle tp is determined from the above, z(s,t) is set by an integration of 
z s = sin(^), where the constant of integration is chosen such that the center of volume of the surface 
x is at the origin. The vertical position z${t) has no bearing on the vesicle shape determination, 
and we hold off further discussion on its dynamics until the following section. 

For the spontaneous curvature model, constraints on the unknown integration length L, the 
surface area A, the volume V, and the two constant Lagrange multipliers P and S are imposed as 

A s = 2vrr, V s = itr 2 sin^, P s = 0, (9) 
£ s = 0, L s = 0. (10) 

Defining Rq as the radius of the sphere with surface area A, the boundary conditions for the five 
constraint equations above are A(0) = V(0) = 0, A(L) = AttRq, V(L) = A-kRqv/3, where v is 
a dimensionless "reduced volume." Due to the Lagrange function being independent of the arc- 
length s, the "Hamiltonian" is a conserved quantity and we have 7(0) = (see Refs. [43, 44]). Also 
defining a reduced spontaneous curvature cq = CqRq, we finally obtain the vesicle morphology as 
set by the two parameters (v, cq). 

In the bilayer coupling model, Eqs. (5-10) are solved with two additional constraints. First, the 
integrated mean curvature M is controlled, M s = n(rK + sin^), and second, a new Lagrangian 
constraint enters, Q s = 0. The system is now closed with boundary conditions on the integrated 
mean curvature: M(0) = and M(L) = 4irRoAa, where Aa is the reduced surface area difference 
between monolayers, Aa = AA/8irRoh. In this case the vesicle morphology is set by the two 
parameters (v, Aa), and the reduced spontaneous curvature c$ has been removed from the shape 
equations via the Legendre transform given above. 

Equations (5-10) are solved numerically. Due to coordinate singularities in the derivatives of 
r and z at the poles, the shape is determined on the contracted interval s 6 [L5,L(1 — 5)] for 
(LS) 1, and Taylor-expanded versions of the boundary conditions are applied. For example, 

r(L5,t) = r(0,t) + (L5)r s (0,t) + O((L5) 2 ) = (L 5) r s (L 5,t) + O ((L5) 2 ) 

= (L 5) cos(V>(L 6)) + O {(L5) 2 ) « Lb. (11) 

To compute the shapes using either model, the arc-length is discretized using m uniformly spaced 
grid points, Sj, with s\ = Lb and s m = L(l — b). A collocation method is then applied in a 
formulation and implementation similar to that recently used by Jiang et al. [44]. We employ a 
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standard continuation scheme in order to interpolate solutions from one point in the parameter 
space (v, co) or (v,Aa) to neighboring points. 

By extremizing the enthalpies F or G, the shape equations give only stationary solutions, not 
necessarily the lowest energy solutions. A numerically determined shape may correspond to an 
energy saddle point, maximum, or minimum. Although it is possible that the lowest energy state 
may not be achievable for a non-equilibrium shape change, for our purposes we will examine the 
minimum energy shapes, and thus a "phase diagram" for the possible shapes is of great use. Just as 
in a more conventional phase transition, shape transformations correspond to transitions between 
different symmetry states. Since we consider only axisymmetric shapes here, spherical solutions 
have the highest symmetry state. For small perturbations around spherical shapes, the solution 
can be represented as 



where the functions YP are the spherical harmonics, and the constants B($ can generate symmetry 
breaking. Because we consider only axisymmetric vesicles, only the m = spherical harmonics 
(of the Y™) contribute to the sum, and the angle 6 is given by tan# = r/z. While it is not 
possible to produce an analytical solution using this formulation, it is useful for understanding the 
morphological transitions in terms of symmetry breaking. For example, breaking t = 2 symmetry 
(-^20 7^ 0) leads to a prolate or oblate shape, while breaking I > 2 symmetry can give more 
complicated shapes, such as the so-called "pear" or "stomatocyte" shapes [7]. In our numerical 
investigation, symmetry is frequently exploited in order to efficiently compute the equilibrium 
shape. In regions of multiple stability the solution branches that correspond to lowest energy 
shapes must be chosen, and by inserting numerically an initial symmetry breaking the algorithm 
used can more readily converge upon the appropriate solution. 

C. Fluid-body interaction 

Modulation of the dimensionless parameter set (v, cq) or (v, Aa) generates quasi-static defor- 
mations which in turn lead to motion in the surrounding fluid medium. Given that the Reynolds 
number is small, the dynamics of the fluid surrounding the vesicle is effectively governed by viscous 
dissipation and is well modeled by the incompressible Stokes equations, 




(12) 



V-cr = 0, V-u = 



(13) 
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where <x = —pi + 2/xE is the Newtonian stress tensor with p the pressure, u the fluid velocity, 
and E the symmetric rate-of-strain tensor, E = |(Vu+ (Vu) T ). The fluid equations are made 
dimensionless by scaling velocities upon U c , lengths upon L c , and time upon t re [ — L c /U c . Since 
the surface area A = AitRq is constant, we define the characteristic length scale by this radius, i.e. 
L c = Rq. Henceforth, the swimming velocity is understood to be dimensionless, and each shape 
cycle occurs over a unit in dimensionless time. 

A no-slip condition is applied on the body surface. For a given path through the parameter 
space (v,co) or (v, Aa), the resulting sequence of instantaneously determined shapes set uniquely 
the "surface deformation velocity" Ud(x, t); namely, 

u d (x(s,0,t),t) = -^(s,<f>,t). (14) 

In addition, the surface moves as a rigid body along the z direction due to axisymmetry, with 
velocity U = Uz = z' (t)z. The no-slip condition is thus written as u(x, t) = Uz, + Ud(x, t). 

To close the system of equations describing the fluid-body interaction, we assume that no 
external forces are acting upon the vesicle, and thus force and torque balance give 

/ tr(x) • n(x) dS = 0, I x x [er(x) • n(x)] dS = 0. (15) 

JS(t) JS(t) 

The computation of the swimming velocity is performed using a standard double-layer boundary 
integral formulation of the Stokes equations. The details of this formulation and numerical method 
are presented in the appendix. 

In addition to computing the swimming velocity, we consider a possibly more important quan- 
tity, the hydrodynamic efficiency. This swimming efficiency is defined as (see Ref. [45]) 



U-F) (U-F 

vh = — f * L = — H ^— r. ( 16 ) 



/ (U + Ud)-fds\ ( / u d -fdS 

JS(t) 1 x Js(t) 



IS® ' x JS(t) 

where f = —a ■ n is the force density acting on the fluid at the body surface, (•) denotes a time- 
average over a full shape cycle, and F = 67r/i a Uz is the force required to move a sphere of radius 
a at a speed U. At each time we use the maximum vesicle radius, a(t) = \\r(s, t)||oo- The first term 
in the denominator of Eq. (16) integrates to zero due to the zero-net force condition (Eq. (15)). 
The computation of the fluid stress a is significantly more involved than the computation of the 
swimming velocity. We employ a numerical method for computing a based on the evaluation 
of a hypersingular integral which may be derived from the double-layer formulation of the fluid 
velocity. The framework and numerical approach are described in the appendix, and a more 
detailed description of the method and examples of its use will be featured in a subsequent paper. 
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Physically, r\n measures the proportion of work done by the vesicle against the surrounding fluid 
which is used for swimming purposes, and is typically on the order of 1% for biological cells. Note 
that the swimming efficiency only measures the hydrodynamic efficiency, not a total efficiency. 
For example, the bending energy of the vesicle is not captured in this measure. The inclusion 
of bending costs into swimming efficiency measures has recently been proposed to study optimal 
locomotion strategies in flagellated cells, but presents an avenue of inquiry beyond the scope of 
this paper [46]. 

IV. VESICLE LOCOMOTION BY SHAPE-CHANGE 

As stated in the introduction, due to the linearity and time-reversibility of Eqs. (13), any 
time-reversible geometrical surface deformations cannot result in a net locomotion. This result is 
known as the Scallop theorem, in reference to the sole, time-reversible motions available to a small 
scallop (opening and closing) [1]. As a consequence of this constraint, a single degree of freedom is 
insufficient for swimming. Two degrees of freedom are however sufficient to generate a swimming 
motion, as first described in Ref. [1], and as we shall show presently for the systems of interest. 

A. Spontaneous curvature model 

We begin by presenting a characteristic shape cycle that can be generated by adjusting the 
reduced volume and spontaneous curvature, (v,co), in a periodic fashion. By selecting a specific 
elliptical path in the (v,co) parameter space, namely v(t) = 0.425 + 0.125 cos(2-7ri), co(t) = —0.1 + 
0.3sin(27rt), the resulting shape cycle is not time-reversible; hence, the constraints of the Scallop 
theorem are bypassed, and locomotion may be achieved. For these parameters the vesicle shapes 
are always stomatocytes, and the neck separating the internal sphere of fluid from the external fluid 
is very small. Figure 3 shows the corresponding minimal energy vesicle shapes at four times, along 
with the vorticity generated in the surrounding fluid by the body deformation, u = V x u. Positive 
vorticity, corresponding to counter-clockwise rotation, is shown in red, and negative vorticity, 
corresponding to clockwise rotation, is shown in blue. Hollow arrows indicate the instantaneous 
swimming velocity of the vesicle, while the plain arrows indicate the direction of time. At zero 
Reynolds number the swimming velocity, external flow, and swimming efficiency are determined 
uniquely by the time-dependent surface geometry and surface deformation velocity, so we need not 
consider the internal flow dynamics (which may in general depend upon the means of modulating 
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FIG. 3: (color online) Stomatocytc shapes and vorticity profiles produced using the spontaneous curvature 
model, with v(t) = 0.425 + 0.125 cos(27rt), Co(t) = —0.1 + 0.3 sin(27rf). Positive vorticity, corresponding to 
counter-clockwise rotation, is shown in red, and negative vorticity, corresponding to clockwise rotation, is 
shown in blue. Hollow arrows indicate the instantaneous swimming velocity. During one cycle, the vesicle 
experience net locomotion in the — z direction. 

the parameters (v, cq)). 

From t = to t = 1/4 the vesicle volume is decreasing while the spontaneous curvature is 
increasing. The decrease in volume draws fluid into the stomatocyte cavity, while the surface 
material near the opening to the cavity moves inward nearly tangentially to the surface itself. 
While the deformation velocity is normal to the surface near the north and south poles (s = and 
s = tt), the deformations are elsewhere primarily tangential, and vorticity is created as the fluid 
is sheared accordingly. At t = 1/2 the vesicle volume is minimal, and the fluid volume inside the 
stomatocyte cavity is beginning to decrease. From t = 1/2 to t = 3/4, the vesicle volume increases 
while the spontaneous curvature continues to decrease to its minimum value. This can best be 
understood by observing that when cq < the membrane prefers a total negative curvature, and 
as can be seen at t = 3/4, the internal cavity of the vesicle takes its smallest value, maximizing 
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FIG. 4: (color online) (a) Phase diagram for the spontaneous curvature model in the (v, cq) parameter space. 
Solid lines are our numerically-calculated lines that denote morphological transitions, while the dashed lines 
are qualitative, and adapted from Rcf. [43]. (b) Three velocity profiles, corresponding to the elliptic paths 
through parameter space indicated in (a), with the largest velocities achieved along the elliptic path enclosing 
the greatest area. 

negative curvature. The increasing volume expels fluid from the cavity, and leads to a reversing of 
the sign of the vorticity. The overall sequence of asymmetric shapes is not time-reversible, leading 
to a net swimming velocity taking place in the — z direction. 

A phase diagram for the minimal energy shapes using the spontaneous curvature model is 
presented in Fig. 4a. The limit lines correspond to discontinuous morphological transitions, and 
therefore cannot be crossed in our quasi-static shape-change approach. One critical line corre- 
sponds to vesicles whose north and south poles self- intersect, and a second line corresponds to 
stomatocyte shapes that have a vanishing opening between the external fluid and the cavity within 
(i.e. the shapes are two spheres, one contained entirely within the other). A third line marks the 
discontinuous phase transition between stomatocyte and oblate shapes. More details may be found 
in Ref. [43]. 

Beyond the symmetry constraints imposed by the Scallop theorem, other symmetry breaking 
is necessary in order for a body to achieve a net motion from a periodic shape cycle. Namely, 
the body surface must express fore-aft asymmetry in order to swim preferentially in any direction. 
Hence, parameter paths in the regions of phase space corresponding to prolate or oblate vesicle 
shapes cannot yield a net motion. However, paths which correspond to stomatocyte or pear shapes 
are fore/aft asymmetric and can swim. Since the area in phase space that contains pear shapes 
is very small, we will only examine the swimming stomatocytes. The largest elliptic path shown 
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in Fig. 4a corresponds to the shape cycle shown in Fig 3. The associated time-evolution of the 
vesicle center of mass velocity is shown in Fig. 4b, along with two other velocities corresponding 
to elliptic paths enclosing smaller areas in Fig. 4a. 

We see in Fig. 4 that the larger the area of the cycle in parameter space, the faster the vesicle 
swims. In fact, the mean velocity roughly scales as the square-root of the area enclosed by the 
elliptic path of phase space. Drawing on an analogy with thermodynamics, cycles with larger area 
in the appropriate ensemble space do more work, and thus we might expect that the transduction of 
shape deformation into mechanical work would exhibit similar behavior. Although our equivalent 
to an equation of state is too complicated to show a simple relationship between swimming velocity 
and the area enclosed in this phase space, the basic idea appears to remain valid. 

We finally note that the net translation during each shape cycle in each case is small compared to 
the amplitude of the motion, and even smaller when compared to the maximum vesicle radius. The 
swimming velocities and hydrodynamic efficiencies of shape cycles in the spontaneous curvature 
model are also small. The maximum velocity achieved for the cycles shown is (U) = —0.008, while 
we calculate an efficiency of r\n = 0.4%. 

B. Bilayer coupling model 

We now consider the bilayer coupling model, for which a schematic phase diagram is shown in 
Fig. 5a. Although in the spontaneous curvature model there are no continuous transitions between 
oblate and stomatocyte shapes, the interesting feature of the bilayer coupling model is the presence 
of a continuous stomatocyte-oblate transition. The upper (solid) line in Fig. 5a denotes a limit line 
between oblate and prolate shapes, while the lower (dashed) line represents a continuous transition 
between stomatocyte and oblate shapes. 

In order to examine how breaking or restoring oblate (£ = 2) symmetry relates to swimming, we 
now consider two shape cycles with equal enclosed area in phase space, as shown in Fig. 5a. The 
upper cycle crosses the continuous transition line, while the lower cycle remains in the stomatocyte 
region. 

The vesicle shapes in the lower cycle of Fig. 5a are displayed in Fig. 6. They correspond to 
a modulation of the volume and surface area difference between monolayers for the vesicle as 
v(t) = 0.775 + 0.075 sin(2vrt), Aa(t) = -0.14 cos(2vrt) + 0.86. From t = to t = 1/4 the vesicle 
volume is increasing, expelling fluid from the cavity and pushing fluid away from the surface of 
the membrane. Due to the larger amount of surface area facing the aft end of the vesicle, the net 
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FIG. 5: (color online) (a) Phase diagram for the bilayer coupling model in the (v, Aa) parameter space, 
adapted from Ref . [43] . The dashed line indicates a continuous transition, while the solid line indicates a limit 
shape. Both lines are shown schematically in order to exaggerate the difference between the shape cycles. 
The two elliptical cycles considered enclose the same area in phase space, but one crosses the transition 
line, (b) Swimming velocity of the vesicle as a function of time, for the two shape cycles shown in (a). 
The squares denote the continuously varying velocity of the lower cycle in (a), which is similar to what we 
observed for the spontaneous curvature model. The circles correspond to the upper cycle in (a) and involves 
a shape transition, and there is a portion of the cycle during which the vesicle has zero swimming velocity 
due to fore-aft symmetry. 



motion during this quarter-cycle is forward. From t = 1/4 to t = 1/2, the "lobes" of the vesicle 
move downwards, propelling the vesicle upwards, albeit at a decreasing rate. This portion of the 
motion resembles the characteristic undulatory shape of a jellyfish, albeit one at zero Reynolds 
number. Between t = 1/2 and t = 3/4, the vesicle deflates and the lobes begin to move upwards 
again, with the material points of the lobes moving almost completely tangentially to the surface. 
This creates a vortex dipole at the lobes, leading to the stagnation point that can be seen in the 
figure. Finally, in the last quarter cycle, the vesicle encloses itself and returns to the starting 
position. We calculate a mean swimming velocity of (U) = —0.048, and a hydrodynamic efficiency 
Of T] H = 0.6%. 

The upper elliptical cycle of Fig. 5a, with shapes illustrated in Fig. 7, follows the parameter 
path v(t) = 0.775 + 0.075 sin(27rf), Aa(t) = -0.14cos(27rt) + 0.89, which lies above the continuous 
stomatocyte-oblate phase transition line from t ~ 0.45 to t ~ 0.55. During this portion of the 
cycle the vesicle has exactly zero swimming velocity due to the fore/aft symmetry of oblate shapes. 
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FIG. 6: (color online) Vesicle shape cycle using the bilayer coupling model, with v(t) = 0.775+0.075 sin(27ri), 
Aa(t) = — 0.14cos(27rf) + 0.89, corresponding to the lower cycle of Fig. 5a. This vesicle does not change 
morphological symmetry states during the swimming cycle and remains within the stomatocyte domain. 
Hollow arrows denote the instantaneous swimming velocity. 

Between t = and t = 1/4, the volume and area difference are decreasing, leading the nearly 
oblate shape into a clearly stomatocyte configuration. For our purposes, we will not address 
the spontaneous symmetry breaking that is associated with crossing a transition line, but simply 
assume that once broken, the cycle will break the symmetry in the same way during each cycle. 
In the example shown, the stomatocyte inflates as it assumes a more oblate shape, expelling fluid 
from the cavity and producing vorticity along the lobes. As the vesicle continues to deflate from 
t = 1/4 to t = 1/2, the lobes sweep downwards, moving the stomatocyte upwards as it assumes 
a perfectly oblate shape. At t ~ 0.42 the shape transitions into an oblate shape, precluding any 
net swimming by symmetry. The swimming velocity as a function of time is shown in Fig. 5b. As 
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FIG. 7: (color online) Vesicle shape cycle using the bilayer coupling model across a continuous phase 
transition, with v(t) = 0.775 — 0.075 sin(27rt), Aa(t) = 0.14 cos(2-7rf) + 0.89, and corresponding to the upper 
cycle of Fig. 5a. This vesicle is oblate for a small part of the cycle, precluding swimming by symmetry, but 
a net locomotion occurs over the entire cycle. Hollow arrows denote instantaneous swimming velocity. 

the oblate vesicle deflates, at t « 0.58 the stomatocyte symmetry state is entered once more, the 
lobes sweep upwards, and the vesicle moves downwards. Despite the presence of a becalmed period 
during the vesicle does not move, the cycle that involves the shape transition yields a larger mean 
velocity than the lower cycle, (U) = —0.055, and an increased hydrodynamic efficiency, t]h = 0.7%. 

As previously noted, crossing the shape transition line between stomatocyte and oblate shapes 
indicated in Fig. 5a yields a continuous shape change. However, if we exploit the analogy with phase 
transitions, we note that some quantities must be discontinuous across the transition. Without 
exploring the details of a dynamic phase transition in the context of vesicle locomotion, although 
the order parameter is continuous, derivatives of the order parameter need not be so. In other 



18 



words, the material at a given point s along the boundary experiences a continuous positional 
change, but a discontinuous velocity relative to the center of mass of the body as the parameters 
are varied continuously through the transition line. The discontinuous relative material velocity 
then generates the discontinuous swimming velocity seen in Fig. 5b for the body which exhibits 
the oblate shapes for part of its periodic cycle. 

Interestingly, even though the area enclosed in phase space by the two cycles illustrated in 
Fig. 5a is the same, the relationship between parameter space, efficiency, and swimming velocity 
is not evident. The upper cycle shown in Fig. 7 has a larger mean swimming speed and is more 
efficient than the cycle shown Fig. 6, suggesting that the vesicle can increase its efficiency by passing 
through a phase transition. 

V. DISCUSSION 

In this paper, we have shown computationally that it is possible for a bilayer vesicle to swim 
under a prescribed shape change using two different vesicle models. By modulating the vesicle 
volume and either its preferred curvature (spontaneous curvature model) or the surface area differ- 
ence between membrane monolayers (bilayer coupling model), the vesicle can be made to undergo 
deformations which are not time-reversible, yielding therefore a net swimming motion. Net lo- 
comotion can be obtained either by continuously modulating fore-aft asymmetric vesicle shapes 
(stomatocytes), or by crossing a continuous shape-transition region with fore-aft symmetric shapes, 
and alternating therefore between fore-aft asymmetric and fore-aft symmetric shapes. 

At first sight, the swimming efficiencies obtained in this paper appear to be low. For the swim- 
ming stomatocyte shown in Fig. 4, the efficiency is on the order of 0.4%, while for the bilayer 
coupling model we calculate an efficiency of 0.6% for a non-transitioning vesicle, and 0.7% for a 
vesicle that undergoes a transition from stomatocyte to oblate. However, it is known from many 
theoretical studies that the hydrodynamic efficiency of swimming microorganisms, such as flagel- 
lated bacteria or spermatozoa, is on the order of 1 to 2% (see Ref. [6] and references therein). Our 
results indicated therefore that the equilibrium morphologies of bilayer vesicles, together with their 
appropriate modulations as is done in this paper, lead to locomotion means which are almost as 
efficient as those displayed by biological cells, and might therefore provide an interesting alternative 
to flagella-based synthetic micro-swimmers. Further optimization of the size and shape of cycle 
in parameter space will likely lead to swimming vesicle outperforming the efficiency of flagellated 
cells. In addition, a swimming vesicle has the advantage that the swimmer and the cargo can be 
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one and the same. 

Let us now discuss the typical time and velocity scales obtained in our simulations. A typical 
vesicle size is approximately 10 /im, and for liposomes k ~ 10~ 19 Nm. Except for very curved 
vesicles, the typical radius of curvature ro is approximately 10 /jm as well, leading to a velocity 
scale of 10 fim/s. This gives calculated mean velocities on the order of 0.1 fim/s for the sponta- 
neous curvature model, and 0.5 /im/s for the bilayer coupling model. Translational and rotational 
diffusion constants for vesicles this size at room temperature are D fa 10~ 14 m 2 /s and D r « 10~ 3 
s _1 , respectively. This implies a time scale for translational diffusion of approximately 10 4 s, and 
a time scale for diffusive reorientation of approximately 10 3 s. Since the actuation proposed in this 
paper can be implemented faster than both of these time scales, significant diffusion will take place 
only after many actuation cycles. For time scales much larger than D^ 1 , the effective vesicle diffu- 
sion will then be given by -D e // ~ U 2 /D r [47], which accounts for both swimming and orientation 
loss. The ratio D e fj/D ~ 10 3 is large, which implies that locomotion will lead to a substantially 
enhanced diffusion of the vesicles over long time scales. 

We have considered only two minimal models for vesicle shape change, and many possible av- 
enues exist to expand upon this basic model, including a study non-axisymmetric vesicles, more 
advanced curvature models, and arc length-dependent spontaneous curvature. Since we have as- 
sumed a quasi-static deformation, non-equilibrium effects would also have to be taken into account 
for fast deformations, and the shape should be fully determined as a balance between elastic and 
fluid forces. In addition, swimming is just one example of behavior that could be exhibited by a 
membrane that is actively deformed. It is perhaps the simplest transduction of geometrical defor- 
mation into mechanical work, and one that we hope provides further inspiration for the combined 
study of membrane physics and low Reynolds number fluid mechanics. 
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Appendix: Velocity and Stress Computation 

The swimming velocity is computed at each time by solving a standard boundary integral 
formulation of the Stokes equations. As an application of the Lorentz reciprocal identity, the 
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solution to Eqs. (13) may be written as integrations upon the surface velocity and the fluid stress, 

u ( x ) = ^7 / G(x, y) • (<r(y) ■ n(y)) dS y + -jj- f u(y) ■ T(x, y) • n(y) dS y (17) 

JS(t) 87T Js(t) 



where 



G -(x,y) = + ^ |^ )(Xj |3 Vl \ (18) 

|x — y| |x — y| d 

T,,(x,y) = _ 6 (^-yO(g-U|)(x fc -I/ fc ) > (19) 

are the singular Stokeslet and Stresslet tensors, respectively (see Ref. [48]). By introducing a 
complementary flow u' which has the same values of the surface force cr ■ n as the flow u on the 
surface S(t), Eq. (17) may be written solely in terms of the second, double-layer integral, 

u(x) = / q(y) • T(x, y) • n(y) dS y , (20) 
JS(t) 

where q(x) is an unknown density of the singular Stresslet tensor. In the limit as the x approaches 
the body surface S(t), inserting the no-slip condition for the surface velocity there, we find the 
expression 

U + u d (x)=/ (q(y)-q(x)).T(x,y).fi(y)dSi,. (21) 
Js(t) 

The vertical swimming velocity U = U • z is related to the Stresslet density as 

4tt 



u A 



[ z ■ q(x) dS (22) 
JS(t) 



IS(t) 

(recall that A is the vesicle surface area). Equation (21) is a well-posed Fredholm integral equation 
of the second kind for the unknown density q(x), and has a unique solution. This approach is 
numerically better conditioned than those based on first-kind equations. 

The Stresslet integral operator in Eq. (21) has a six-dimensional nullspace corresponding to rigid 
body motion, and in the presence of external body forces or torques this representation must be 
closed by a range completion technique (see Ref. [49]). However, in the swimming problem where 
the deformation velocity Ud(x) is specified and there are no body forces or torques, Eqs. (21-22) 
are closed and prescribe uniquely the swimming velocity U. 

The integrand in Eq. (21) is discontinuous at the singularity but finite, so that the integrals are 
computed to second-order in the surface mesh element size using a standard trapezoidal quadrature 
(setting the quadrature weight to zero at the singularity). The axisymmetry of the problem is 
inserted into the definition of the body surface as well as the density q(x) . The number of gridpoints 
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is chosen to be sufficiently large such that further resolution does not significantly alter the density 
q(x) or the swimming velocity U. 

At each time, the curve (r(s,t), z(s,t)) is discretized uniformly in s. Application of a Nystrom 
collocation method produces a linear system of equations for the density q(x) at the gridpoints, 
which is then solved iteratively using the method GMRES [50], with an inversion error tolerance 
such that the only errors are due to discritization. Finally, the body position Zo(t) is updated at 
each time using a second-order Runge-Kutta method. Both convergence tests and comparison with 
known exact solutions were used to validate the code [51-54]. 

Computing the hydrodynamic or swimming efficiency (which requires pointwise information 
about the stress <r) is more difficult. Here we compute <x(x) using the approach outlined below, 
though a more detailed description of the method and examples of its use will be featured in a 
subsequent paper. 

Many common methods for computing the stress are developed using a first-kind boundary 
integral formulation of the Stokes equations, and hence can suffer from the ill-posedness of the 
underlying equations [48]. Instead, we solve for the surface stress by evaluating a hypersingular 
integral which may be derived from the second- kind integral equation for the velocity (see Ref. [48]), 



/' 

where 



-cj im (x) = / qj(y)£ ijkm (x,y)nk(y)dSy, (23) 
W Js 



r ( \- A dim^jk R {xk - yk)[5jm(Xj ~ Pi) + 5jj(x T 
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|x-y| J |x-yT 

_ g pEj ~ yj)[5mk{Xi - Vi) + SjkjXrn ~ Vrn)] QQ (xj ~ Vi){Xj ~ Vj)(x k - yk){x m ~ Vm) 

|x-y| 5 |x-y| 7 

and we have set S{t) = S for clarity. The expression £(x, y) is achieved by differentiating the 
double-layer integral for the fluid velocity and including the pressure term which may also be 
written as an integration against q(x), with a = —pi + /i (Vu + Vu r ) (see Ref. [48]). The stress 
is determined on the same spatial grid as used to determine the swimming velocity, a uniform 
discritization in s of the curve (r(s, t), z(s, t)) (with polar angle (ft = 0). The integration of Eqn (23) 
is performed in local polar coordinates, and the singular contributions are handled analytically as 
follows. The procedure follows the work of Guiggiani et al. [55]. 

The integration of Eq. (23) is performed on a modified surface S = s t + (S — e e ) and is taken in 
two parts: the portion of a sphere of radius e centered at the singular point x which is internal to 
the body surface (s e ) and intersects the surface S at its boundary, and the body surface punctured 
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by the sphere (S — e e ). The modified surface limits to the body surface S as e — > 0. For a point 
x£S, Eq. (23) is written ELS £L small e limit, 



lm il / qj(y)£ijkm(x,y)nk(y)dS y + q j (y)C ijkm (x,y)h k (y)dS y \- (25) 



Under the assumption that q(x) is differ entiable, with a derivative which is Holder continuous, we 
subtract and add the density q(x) and its gradient at the singular point in the second integral of 
Eq. (25), 

-o-im(x) = lim{ / q j (y)C ijkm (x,y)h k (y)dS y (26) 

+ (?j(y) - ?j( x ) - (zfc - yh)qj,h(x)^£ijkm(x,y)n k (y) ds y (27) 
+ ?i,ft( x ) / (^h - yft)Ajfem( x , y)^fc(y) ^ (28) 

+ qj(x) C ijkm (x,y)n k (y)dS y ,\ (29) 

where q^/j = dqj/dxh- As shown in Ref. [55], the above integration may be reduced to a final 
formula upon the introduction of a local polar coordinate system (p, rj) about the target point 
x(s, (/>), with 

4> = (j) + P cos (77), s' = s + psin(?/), (30) 

where 77 € [0, 27r), p € [0, p(r/)], and 

dS'y = J(s')ds' dtp' = J(s'(p,r]))pdpdr], (31) 

with J(s') = |x^ x x^| the surface Jacobian. p = p(rj) is the equation in the local polar coordinate 
system of the edge of the semi-periodic domain, (s,<fi) G ([0, L] x [0, 2tt]). The integration is assisted 
by the extra factor of p in the surface area element, and the final expression for the fluid stress 
may be reduced to 



a im (x) = j J { Fijk (p, V )-[-^^ + ^^]}dpdr, (32) 
{4^ fo) ^ - F^fr) [-L] } ^, (33) 



/' 

where Fy fc (x,y) = %(x)£yfc(x,y)n fe (y) [55]. The functions F^ 1 ' '(77) and F^~ 2 \r]) are the singular 
parts of an expansion of Fyfc(p, 77) about p = 0. The integrals above all have finite integrands, and 
are treated using adaptive quadrature methods. 
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Convergence tests and comparisons with known exact solutions were used to validate the code. 
In particular, we have checked to ensure that the surface deformation relation of Samuel & Stone 
(1996) is satisfied [52, 54]. With the stress <x in hand, the efficiency tjh (Eq. (16)) is determined 
to second-order in the grid-spacing by a simple trapezoidal quadrature. The stress need only be 
computed for = due to axisymmetry. 

As a final note, at zero Reynolds number the swimming velocity and efficiency are entirely 
determined by the surface deformation velocity. Other more general measures of energetic expen- 
diture and total efficiency have been considered for other swimming systems (see Ref. [46]), but in 
this case the total efficiency will depend significantly upon the means used to produce the vesicle 
shape-change. In addition, should there be a fluid internal to the vesicle, for example, internal 
dissipation costs would be relevant in a more general measure of energetic expenditure. 
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